Using only tracks that passed QC
res = lapply(1:4, function(d) {
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:d), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
p$Day = d
return(p)
})
res = do.call('rbind', res)
res = left_join(res, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])
# Rank day 4
p2 = res %>% filter(Day == 4) %>% group_by(GroupNumber, Sex) %>%
mutate(
Rank = rank(-DS_cum),
Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
)
res = left_join(res, p2[,c('RealMouseNumber', 'Rank')])r2 = res %>% group_by(Sex, Day, Rank) %>%
dplyr::summarize(
M = mean(DS_cum, na.rm = T),
se = sd(DS_cum, na.rm = T)/sqrt(n()),
)
f1_1 = ggplot(r2, aes(x = as.factor(Day), y = M, fill = Rank,
color = Rank)) +
geom_line(data = res, aes(group = RealMouseNumber, y = DS_cum), alpha = 0.4, size = 0.4) +
geom_errorbar(aes(ymin = M-se, ymax = M+se), width = 0) +
geom_line(aes(group = Rank), size = 0.8) +
labs(y = "Baseline David's Score", x = 'Day') +
geom_point(size = 2.5, pch = 21, color = 'grey20') +
scale_fill_manual(values = obj$colors_rank, name = 'Final rank',
labels = c('\u03B1', '\u03B2', '\u03B3', '\u03B4')) +
scale_color_manual(values = obj$colors_rank, name = 'Final rank',
labels = c('\u03B1', '\u03B2', '\u03B3', '\u03B4')) +
scale_x_discrete(expand = c(0.1, 0.1)) +
facet_grid(~Sex) + theme(legend.position = 'top')
f1_1wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = wins, nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F)
h = h %>% dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID)))
h = left_join(h, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])
p = subset(obj$profiles, !Suspect & Day %in% 2:4) %>%
dplyr::select(RealMouseNumber, all_props) %>%
group_by(RealMouseNumber) %>%
dplyr::summarise_at(.vars = all_props, mean, na.rm = T)
d = left_join(h, p)
idVars = c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props = c('DistanceOutside', 'GridEntropy6x6', 'FractionOfTimeOutside')
p <- d %>%
dplyr::select(idVars, props) %>%
gather(prop, value, -idVars) %>%
mutate(prop = plyr::mapvalues(prop, all_props, all_props_pretty))f1_2 = ggplot(p, aes(x = DS, y = value, fill = Sex)) +
geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
geom_point(pch = 21, size = 2) +
labs(x = "Baseline David's Score", y = '') +
facet_grid(prop~Sex, scales = 'free_y') +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
theme(legend.position = 'none')
f1_2wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = wins, nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F)
h = h %>% dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID)))
h = left_join(h, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])
p = subset(obj$profiles, !Suspect & Day %in% 2:4) %>%
dplyr::select(RealMouseNumber, all_props) %>%
group_by(RealMouseNumber) %>%
dplyr::summarise_at(.vars = all_props, mean, na.rm = T)
d = left_join(h, p)
idVars = c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props = all_props
p <- d %>%
dplyr::select(idVars, props) %>%
gather(prop, value, -idVars) %>%
group_by(prop, Sex) %>% mutate(N = n()) %>%
group_by(prop, Sex, N) %>% nest() %>%
mutate(
M = map(data, function(d) {
tidy(cor.test(d$value, d$DS, method = 'spearman'))
})
) %>% dplyr::select(-data) %>% unnest() %>%
mutate(
prop = fct_reorder(prop, abs(estimate))
) %>% group_by(Sex) %>%
mutate(qval = p.adjust(p.value, method = 'BH')) %>%
arrange(p.value) %>%
mutate(sig_qval = qval <= 0.1) Social box features with significant q-values:
| Female | Male | |
|---|---|---|
| MedianContactDuration | 0 | 0 |
| MeanContactDuration | 0 | 0 |
| FractionOfTimeInTheOpenOutside | 0 | 0 |
| MeanSpeedOutside | 0 | 0 |
| FractionOfTimeInSmallNestOutside | 0 | 0 |
| FractionOfTimeInFeederOutside | 0 | 0 |
| FoodOrWaterPerTimeOutside | 0 | 0 |
| MedianSpeedOutside | 0 | 0 |
| FractionOfTimeOnRampOutside | 0 | 0 |
| DistanceFromWallsInOpen | 0 | 0 |
| GridMI6x6 | 0 | 0 |
| EntropyOutside | 0 | 0 |
| DistanceFromNest | 0 | 0 |
| HighPlacePerTimeOutside | 0 | 0 |
| ContactRate | 0 | 0 |
| VisitsOutsideRate | 0 | 0 |
| BeingApproachedRate | 0 | 0 |
| ForagingCorrelation | 0 | 0 |
| AngularVelocity | 0 | 0 |
| FractionOfTimeAtHighPlace | 0 | 0 |
| FractionOfTimeInWaterOutside | 0 | 1 |
| ApproachRateOutside | 1 | 0 |
| TangentialVelocity | 0 | 0 |
| DiffBetweenApproachesAndChases | 0 | 0 |
| FractionOfTimeInContactOutside | 0 | 0 |
| ContactRateOutside | 0 | 1 |
| FractionOfTimeNearFoodOrWater | 0 | 1 |
| DistanceOutside | 0 | 1 |
| FractionOfTimeInLabyrinthOutside | 1 | 0 |
| NumberOfApproachs | 0 | 1 |
| ApproachRate | 0 | 1 |
| FractionOfTimeAloneOutside | 0 | 1 |
| ProximateVsDistantFood | 0 | 1 |
| GridEntropy6x6 | 0 | 1 |
| Entropy | 0 | 1 |
| NumberOfApproachesPerMiceOut | 0 | 1 |
| FractionOfTimeOutside | 0 | 1 |
| BeingApproachedRateOutside | 0 | 1 |
| FractionOfBeingApproachedPerContact | 0 | 1 |
| FractionOfApproachEscapeBehaviorPerAggression | 0 | 1 |
| AggressiveEscapeRate | 0 | 1 |
| FractionOfApproachesPerContact | 0 | 1 |
| BeingFollowedRate | 1 | 1 |
| NAEscapeRate | 1 | 1 |
| AggressiveChaseRateOutside | 1 | 1 |
| AggressiveEscapeRateOutside | 0 | 1 |
| FractionOfEscapesPerContact | 0 | 1 |
| AggressiveChaseRate | 1 | 1 |
| NAEscapeRateOutside | 1 | 1 |
| FractionOfChasesPerContact | 1 | 1 |
| NAChaseRateOutside | 1 | 1 |
| FollowRateOutside | 1 | 1 |
| BeingFollowedRateOutside | 1 | 1 |
| FollowRate | 1 | 1 |
| NAChaseRate | 1 | 1 |
| FractionOfBeingFollowedPerContact | 1 | 1 |
| FractionOfNAEscapesPerContact | 1 | 1 |
| FractionOfFollowsPerContact | 1 | 1 |
| FractionOfNAChasesPerContact | 1 | 1 |
| ProximateVsDistantWater | 0 | 0 |
Distance Outside:
Fraction of Time Outside:
Grid Entropy [6x6]
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>%
dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(MouseID)) %>%
group_by(GroupNumber, Sex, Condition, Day) %>%
mutate(DS_rank = rank(-DS))
p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')
p = p %>% group_by(GroupNumber, Sex, Condition, MouseID) %>% arrange(GroupNumber, Sex, MouseID, Day) %>%
dplyr::select(-DS) %>% rename(DS = DS_rank) %>%
mutate(
lag0 = DS, lag1 = lag(DS), lag2 = lag(DS, 2), lag3 = lag(DS, 3), lag4 = lag(DS, 4),
lead0 = DS, lead1 = lead(DS), lead2 = lead(DS, 2), lead3 = lead(DS, 3), lead4 = lead(DS, 4)
)
pp = p %>% group_by(Day, GroupNumber, Sex) %>% nest() %>%
mutate(
lead = map(data, function(d) {
x = cor(d[, c('lead0', 'lead1', 'lead2', 'lead3')],
use = "pairwise.complete.obs", method = 'pearson') %>%
melt() %>% filter(Var1 == 'lead0')
})
) %>% dplyr::select(-data) %>% unnest() %>%
mutate(uGroup = paste0(GroupNumber, '_', Sex))pd = 0.2
pp %>% group_by(Var2, Sex, Day) %>%
dplyr::summarize(
M = mean(value, na.rm = T),
se = sd(value, na.rm = T)/sqrt(n())) %>%
filter(Day == 1) %>% ungroup() %>%
mutate(Var2 = plyr::mapvalues(Var2, c('lead0', 'lead1', 'lead2', 'lead3'), c(1:4))) %>%
ggplot(aes(x = Var2, y = M, color = Sex, fill = Sex)) +
# geom_hline(yintercept = 0, lty = 2, color = 'grey') +
geom_line(aes(group = Sex), position = position_dodge(pd)) +
ylab("David's Score\ncorrelation with day 1") + xlab('Day') + facet_grid(~Day, scales = 'free_x') +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
geom_errorbar(aes(ymax = M+se, ymin = M-se), width = 0, position = position_dodge(pd)) +
geom_point(size = 2, pch = 21, color = 'grey20', position = position_dodge(pd)) +
theme(legend.position = 'top')# Generate a distribution of expected rank changes
set.seed(42)
sampleGroupRanks = function (i) {sample(1:i, size = i, replace = F)}
# Fake data
nMice = length(unique(obj$master$RealMouseNumber))
nDays = 4
nIter = 10
ref = lapply(1:nIter, function(i) {
dat = data.frame(
Day = rep(1:nDays, nMice),
MouseNumber = rep(1:nMice, each = nDays),
GroupNumber = rep(1:(nMice/4), each = 4*nDays)
)
dat = dat %>% group_by(GroupNumber, Day) %>%
nest() %>%
mutate(
Rank = map(data, function(d) sampleGroupRanks(nrow(d)))
) %>% unnest()
d = dat %>% group_by(GroupNumber, MouseNumber) %>%
arrange(GroupNumber, MouseNumber, Day) %>%
mutate(
lagRank = lag(Rank, 1),
RankChange = abs(Rank - lagRank)
) %>% filter(!is.na(lagRank))
ref = d %>% group_by(GroupNumber, MouseNumber, RankChange) %>%
tally() %>%
group_by(GroupNumber, MouseNumber) %>%
dplyr::mutate(
s = sum(n),
SumRankChangesFraction = n / sum(n)
) %>%
group_by(RankChange) %>%
tally() %>%
dplyr::mutate(
SumRankChangesFraction = n / sum(n)
)
ref$iter = i
return(ref)
})
ref = do.call('rbind', ref)
ref = ref %>% group_by(RankChange) %>%
dplyr::summarise(
M = mean(SumRankChangesFraction, na.rm = T),
Qt = qt(0.975,n()-1) * sd(SumRankChangesFraction)/sqrt(n()),
CI_95 = M + Qt,
CI_05 = M - Qt
)
# Real final ranks
r = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
r = left_join(r, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])
r = r %>% group_by(GroupNumber, Sex) %>%
mutate(
Rank = rank(-DS_cum),
Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
)
# Real ranks throughout
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>%
dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(MouseID)) %>%
group_by(GroupNumber, Sex, Condition, Day) %>%
mutate(DS_rank = rank(-DS))
p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')
p = left_join(p, r[,c('RealMouseNumber', 'Rank')])
p = p %>% group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>%
arrange(GroupNumber, Sex, MouseID, Day) %>%
mutate(
lagRank = lag(DS_rank),
RankChange = abs(DS_rank - lagRank)
) %>% filter(!is.na(RankChange))
pp = p %>%
# group_by(Sex, Rank, RankChange, RealMouseNumber) %>%
group_by(Sex, Rank, RankChange) %>%
tally() %>%
# group_by(Sex, Rank, RealMouseNumber) %>%
group_by(Sex, Rank) %>%
mutate(
s = sum(n),
SumRankChangesFraction = n / sum(n, na.rm = T)
)
ggplot(pp, aes(x = RankChange, y = SumRankChangesFraction)) +
geom_ribbon(data = ref, aes(ymax = CI_95, ymin = CI_05, y = M),
fill = 'grey', alpha = 0.5) +
geom_line(data = ref, color = 'grey', lty = 1, aes(y = M)) +
geom_line(aes(col = Rank)) +
geom_point(aes(fill = Rank), pch = 21, color = 'grey20', size = 2.5) +
facet_grid(.~Sex) +
labs(y = 'Probability') +
scale_color_manual(values = obj$colors_rank) +
scale_fill_manual(values = obj$colors_rank) +
theme(legend.position = 'top')x = left_join(pp, ref)
b = seq(-2, 2, 1)
l = c(-4, -2, 0, 2, 4)
ggplot(x, aes(x = as.factor(RankChange), y = log2(SumRankChangesFraction / M))) +
geom_hline(yintercept = 0, lty = 2, color = 'grey') +
geom_line(aes(col = Rank, group = Rank), size = 0.5,
position = position_dodge(0.4)) +
geom_point(aes(fill = Rank, col = Rank), pch = 21,
color = 'grey20',
size = 2.5,
position = position_dodge(0.4)) +
labs(y = 'Fold change in probability\n(true / expected)', x = 'Absolute rank change') +
facet_grid(.~Sex) +
scale_y_continuous(breaks = b, labels = l) +
scale_fill_manual(values = obj$colors_rank) +
scale_color_manual(values = obj$colors_rank) +
theme(legend.position = 'none')# Real final ranks
r = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
r = left_join(r, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])
r = r %>% group_by(GroupNumber, Sex) %>%
mutate(
Rank = rank(-DS_cum),
Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
)
# Real ranks throughout
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>%
dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(MouseID)) %>%
group_by(GroupNumber, Sex, Condition, Day) %>%
mutate(DS_rank = rank(-DS))
p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')
p = left_join(p, r[,c('RealMouseNumber', 'Rank')])
p = p %>% group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>%
arrange(GroupNumber, Sex, MouseID, Day) %>%
mutate(
lagRank = lag(DS_rank),
RankChange = abs(DS_rank - lagRank)
) %>% filter(!is.na(RankChange))
pp = p %>%
dplyr::mutate(
NoChange = RankChange == 0
)
ppp = pp %>% group_by(Rank, Sex) %>%
dplyr::summarise(
N = n(),
S = sum(NoChange)
) %>%
dplyr::mutate(
P = S / N
)
ppp$sig = NA
ppp$sig[ppp$Rank == 'alpha'] = '***'
ppp$sig[ppp$Rank == 'delta' & ppp$Sex == 'Male'] = '***'
ppp$Rank = plyr::mapvalues(ppp$Rank, c('alpha', 'beta', 'gamma', 'delta'),
c( '\u03B1', '\u03B2', '\u03B3', '\u03B4'))
f1_3 = ggplot(ppp, aes(x = Rank, y = P / .25, fill = Rank)) +
geom_bar(stat = 'identity', width = 0.7, alpha = 1, aes(color = Rank)) +
geom_hline(yintercept = 1, lty = 1, color = 'grey30') +
geom_text(aes(label = Rank, y = 1), nudge_y = 0.08, parse = F, fontface = 'italic', color = 'grey30') +
geom_text(aes(label = N), color = 'grey20', nudge_y = 0.08, parse = F, size = 3, fontface = 'italic', color = 'grey30') +
geom_text(aes(label = sig), color = 'grey20', nudge_y = 0.15, parse = F, size = 4, fontface = 'bold', color = 'grey30') +
facet_grid(.~Sex) +
labs(y = 'Rank maintenance odds', x = '') +
scale_color_manual(values = obj$colors_rank) +
scale_fill_manual(values = obj$colors_rank) +
scale_y_continuous(trans = 'log2', breaks = c(0.66, 1, 1.5, 2, 2.5, 3), limits = c(0.66, 3)) +
theme(legend.position = 'none',
axis.text.x = element_blank())
f1_3 Stats:
res = lapply(1:nrow(ppp), function(i) {
d = ppp[i,]
tidy(binom.test(d$S, n = d$N, p = .25, alternative = 'greater'))
})
res = do.call('rbind', res)
ppp = cbind(as.data.frame(ppp), as.data.frame(res)) %>% arrange(Sex, Rank)
as.data.frame(ppp)baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest()
stress = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 5), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>%
rename(DS_stress = DS)
h = left_join(baseline, stress)Baseline v. Stress DS
f1_4 = ggplot(h, aes(x = DS, y = DS_stress, fill = Sex)) +
geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
geom_point(pch = 21, size = 2) +
labs(y = "David's Score\nfollowing acute stress", x = "Baseline David's Score") +
facet_grid(.~Sex) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
theme(legend.position = 'none')
f1_4Males:
Females
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 2:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>%
dplyr::select(GroupNumber, Day, Sex, Condition, data) %>% unnest() %>%
group_by(GroupNumber, Sex, Condition, MouseID, Day) %>%
dplyr::summarize(wins = sum(wins, na.rm = T)) %>% # Collapse over OtherID
group_by(GroupNumber, Sex, Condition, MouseID) %>%
dplyr::summarize(wins = mean(wins, na.rm = T)) # Collapse over Day
stress = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 5), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, data) %>% unnest() %>%
rename(wins_stress = wins) %>%
group_by(GroupNumber, Sex, MouseID) %>%
dplyr::summarize(wins_stress = sum(wins_stress, na.rm = T)) # Collapse over OtherID
h2 = left_join(baseline, stress)Chases from baseline to stress
pp = h2 %>% group_by(GroupNumber, MouseID, Sex, Condition) %>%
summarise_at(.vars = c('wins', 'wins_stress'), mean, na.rm = T) %>%
gather(Stress, Chases, -GroupNumber, -Sex, -Condition, -MouseID) %>%
mutate(Stress = plyr::mapvalues(Stress, c('wins', 'wins_stress'), c('Baseline', 'Acute Stress')),
Stress = factor(Stress, levels = c('Baseline', 'Acute Stress')))
Md = pp %>% group_by(Sex, Stress) %>%
dplyr::summarise(
Md = median(Chases, na.rm = T),
Mn = mean(Chases, na.rm = T),
SE = sd(Chases, na.rm = T) / sqrt(n()),
MAD = stats::mad(Chases, na.rm = T)
)
# ppp = left_join(pp, Md)
f1_5 = ggplot(pp,
aes(x = Chases, fill = Sex, color = Sex)) +
geom_density(alpha = 0.7) +
geom_rug(aes(color = Sex), size = 1, show.legend = F) +
geom_errorbarh(data = Md[Md$Sex == 'Male',], inherit.aes = F,
aes(xmin = Md, xmax = Md + MAD, y = -0.1, col = Sex), height = 0.1, show.legend = F) +
geom_point(data = Md[Md$Sex == 'Male',],
pch = 21, size = 2.5, aes(x = Md, y = -0.1), color = 'grey20', show.legend = F) +
geom_errorbarh(data = Md[Md$Sex == 'Female',], inherit.aes = F,
aes(xmin = Md, xmax = Md + MAD, y = -0.2, col = Sex), height = 0.1, show.legend = F) +
geom_point(data = Md[Md$Sex == 'Female',],
pch = 21, size = 2.5, aes(x = Md, y = -0.2), color = 'grey20', show.legend = F) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
facet_grid(Stress~.) +
scale_x_log10() + labs(y = '') +
theme(legend.position = c(0.1, 0.9),
legend.title = element_blank(),
panel.grid.major.x = element_line(color = 'grey', size = 0.2, linetype = 2),
# panel.grid.minor.x = element_line(color = 'grey', size = 0.1),
panel.border = element_blank(),
axis.text.y = element_blank())
f1_5s = h2 %>% gather(Stage, Chases, -GroupNumber, -Sex, -Condition, -MouseID)
s$Stage = plyr::mapvalues(s$Stage, c('wins', 'wins_stress'), c('Baseline', 'Stress'))
s = left_join(s, obj$mice[,c('GroupNumber', 'Sex', 'MouseID', 'RealMouseNumber')])
summary(aov(log(Chases + 1) ~ Sex * Stage + Error(as.factor(RealMouseNumber)), data = s))##
## Error: as.factor(RealMouseNumber)
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.34 0.3366 0.211 0.648
## Stage 1 0.59 0.5950 0.372 0.544
## Sex:Stage 1 0.01 0.0093 0.006 0.940
## Residuals 80 127.88 1.5984
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 5.704 5.704 30.559 0.000000415 ***
## Sex:Stage 1 0.004 0.004 0.023 0.881
## Residuals 78 14.559 0.187
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
col1 = cowplot::plot_grid(plotlist = list(f1_1, f1_2), labels = c('a', 'b'), label_fontface = 'plain',
ncol = 1, align = 'v', axis = 'lr', rel_heights = c(1, 2))
col2 = cowplot::plot_grid(plotlist = list(f1_3, f1_4, f1_5), labels = c('c', 'd', 'e'), label_fontface = 'plain',
ncol = 1, align = 'v', axis = 'lr', rel_heights = c(1, 0.8, 1.2))
cowplot::plot_grid(plotlist = list(col1, col2), ncol = 2, align = 'h')if(exportFigures) {
ggsave('../manuscript/figures/1_Panel.pdf', width = 9.5, height = 9.3, dpi = 1200)
ggsave('../manuscript/figures/1_Panel.png', width = 9.5, height = 9.3, dpi = 1200)
}dat = obj$tests_combined_adj
dat$Sex = factor(dat$Sex, levels = c('Female', 'Male'))
dat$SexCondition = paste0(dat$Sex, '_', dat$Condition)
dat$SexCondition = factor(dat$SexCondition,
levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))Stats:
dat$value = dat[['ORG_BodyweightChange']]
dat %>% group_by(Sex, Condition) %>%
dplyr::summarize(
n = n(),
n_groups = length(unique(GroupNumber)),
M = mean(value, na.rm = T),
sd = sd(value, na.rm = T)
)dat %>% group_by(Sex, Condition) %>% nest() %>%
mutate(SW = map(data, function(d) shapiro.test(d$value)),
tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()Normality is not violated in any group
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 3 | 2.658661 | 0.0550054 |
| 69 | NA | NA |
Groupwise variances are not heteroscedastic (although a trend exists)
Males
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 1 | 3.419868 | 0.0742915 |
| 30 | NA | NA |
Females
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 1 | 2.793313 | 0.1026649 |
| 39 | NA | NA |
Pairwise
dat %>% group_by(Sex) %>% nest() %>%
dplyr::mutate(
t = map(data, function(d) {
x = tidy(t.test(d$value ~ d$Condition))
# x = tidy(kruskal.test(value ~ Condition, d))
}
)
) %>% dplyr::select(-data) %>% unnest()dd = data.frame(Sex = c('Female', 'Male'), Sig = c('***', '*'), y = c(20, 24))
f2_1 = ggplot(dat, aes(x = Condition, y = ORG_BodyweightChange, fill = SexCondition)) +
geom_hline(yintercept = 0, color = 'grey', lty = 2) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
geom_point(pch = 21, size = 2,
position = position_jitterdodge(jitter.width = 0,
jitter.height = 0, dodge.width = 0.5)) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
scale_y_continuous(breaks = seq(-30, 30, 10)) +
labs(x = '', y = 'Bodyweight change') +
facet_grid(.~Sex) + theme(legend.position = 'none') +
geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_1d2 = dat %>% group_by(Sex, Condition, SexCondition) %>%
dplyr::summarize(
Md = median(ORG_FurScore_Cumulative, na.rm = T),
MAD = stats::mad(ORG_FurScore_Cumulative, na.rm = T)
)
dd = data.frame(Sex = c('Female', 'Male'), Sig = c('***', '***'), y = c(4, 6))
f2_2 = ggplot(dat, aes(x = Condition, y = ORG_FurScore_Cumulative, fill = SexCondition)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
labs(x = '', y = 'Coat state') +
facet_grid(.~Sex) + theme(legend.position = 'none') +
geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_2Stats:
dat$value = dat[['ORG_FurScore_Cumulative']]
dat %>% group_by(Sex, Condition) %>%
dplyr::summarize(
n = n(),
n_groups = length(unique(GroupNumber)),
M = mean(value, na.rm = T),
sd = sd(value, na.rm = T)
)dat %>% group_by(Sex, Condition) %>% nest() %>%
mutate(SW = map(data, function(d) shapiro.test(d$value)),
tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()Normality is violated in both male groups as well as female controls
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 3 | 1.535438 | 0.2130927 |
| 69 | NA | NA |
Groupwise variances are not heteroscedastic
dat = dat %>%
ungroup() %>%
# dplyr::mutate(value_trans = sqrt(value + 1)) %>%
dplyr::mutate(value_trans = rank(value))
ggplot(dat, aes(x = Condition, y = value_trans, fill = SexCondition)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
labs(x = '', y = 'Coat state') +
facet_grid(.~Sex) + theme(legend.position = 'none')dat %>% group_by(Sex, Condition) %>% nest() %>%
mutate(SW = map(data, function(d) shapiro.test(d$value_trans)),
tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()No transformation helped for the female controls
Condition:
Sex:
Condition per sex
dat %>% group_by(Sex) %>%
nest() %>%
dplyr::mutate(
t = map(data, function(d) {
tidy(kruskal.test(value ~ Condition, d))
})
) %>% dplyr::select(-data) %>% unnest()dd = data.frame(Sex = c('Female', 'Male'), Sig = c('', '**'), y = c(0, 0))
f2_3 = ggplot(dat, aes(x = Condition, y = ORG_Adrenals_adj, fill = SexCondition)) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
geom_point(pch = 21, size = 2,
position = position_jitterdodge(jitter.width = 0,
jitter.height = 0, dodge.width = 0.5)) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
labs(x = '', y = 'Adrenal weight') +
facet_grid(.~Sex) + theme(legend.position = 'none') +
geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_3Stats:
dat$value = dat[['ORG_Adrenals_adj']]
dat %>% group_by(Sex, Condition) %>%
dplyr::summarize(
n = n(),
n_groups = length(unique(GroupNumber)),
M = mean(value, na.rm = T),
sd = sd(value, na.rm = T)
)dat %>% group_by(Sex, Condition) %>% nest() %>%
mutate(SW = map(data, function(d) shapiro.test(d$value)),
tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()Normality is not violated in any group
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 3 | 0.4508866 | 0.7174833 |
| 69 | NA | NA |
Groupwise variances are not heteroscedastic
Test
Test
Without interaction:
Pairwise
dat %>% group_by(Sex) %>% nest() %>%
dplyr::mutate(
t = map(data, function(d) {
x = tidy(t.test(d$value ~ d$Condition))
# x = tidy(kruskal.test(value ~ Condition, d))
}
)
) %>% dplyr::select(-data) %>% unnest()f2_4 = ggplot(dat, aes(x = Condition, y = ORG_Cort, fill = SexCondition)) +
# geom_hline(yintercept = 0, color = 'grey', lty = 2) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
geom_point(pch = 21, size = 2,
position = position_jitterdodge(jitter.width = 0,
jitter.height = 0, dodge.width = 0.5)) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
labs(x = '', y = 'Corticosterone') +
facet_grid(.~Sex) + theme(legend.position = 'none')
f2_4Stats:
dat$value = dat[['ORG_Cort']]
dat %>% group_by(Sex, Condition) %>%
dplyr::summarize(
n = n(),
n_groups = length(unique(GroupNumber)),
M = mean(value, na.rm = T),
sd = sd(value, na.rm = T)
)dat %>% group_by(Sex, Condition) %>% nest() %>%
mutate(SW = map(data, function(d) shapiro.test(d$value)),
tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()Normality is violated in all groups
Transform the data didn’t help
Levene’s Test
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 3 | 3.555961 | 0.0186511 |
| 69 | NA | NA |
Groupwise variances are heteroscedastic
Per sex Males - not heteroscedastic
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 1 | 0.4056222 | 0.5290303 |
| 30 | NA | NA |
Females - not heteroscedastic
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 1 | 0.2798986 | 0.5997662 |
| 39 | NA | NA |
Condition:
Sex (not fair, since variances are unequal):
Condition per sex
dat %>% group_by(Sex) %>%
nest() %>%
dplyr::mutate(
t = map(data, function(d) {
tidy(kruskal.test(value ~ Condition, d))
})
) %>% dplyr::select(-data) %>% unnest()dat = obj$tests_combined_adj
props = colnames(dat)[regexpr('_', colnames(dat)) > 0]
props_pretty = c(
'EPM | Center time',
'EPM | Closed arm distance',
'EPM | Closed arm entries',
'EPM | Closed arm mean visit',
'EPM | Closed arm time',
'EPM | Distance',
'EPM | Open arm distance from center',
'EPM | Open arm distance',
'EPM | Open arm entries',
'EPM | Open arm mean visit',
'EPM | Open arm time',
'EPM | Closed arm preference',
'NBT | Intact material',
'NBT | Nest score',
'OFT | Mean distance from center',
'OFT | Distance',
'OFT | Inner zone mean speed',
'OFT | Inner zone distance',
'OFT | Inner zone visits',
'OFT | Inner zone mean visit duration',
'OFT | Max. speed',
'OFT | Meander',
'OFT | Outer zone mean speed',
'OFT | Outer zone distance',
'OFT | Outer zone mean visit duration',
'OFT | Outer zone preference',
'OFT | Immobility time',
'PHYS | Adrenal weight',
'PHYS | Bodyweight change',
'PHYS | Corticosterone',
'PHYS | Coat state',
'SPT | Median preference',
'SPT | Total preference',
'SPL | Grooming time',
'SPL | Latency to groom',
'TST | Immobility episodes',
'TST | Immobility time'
)
as.data.frame(cbind(props, props_pretty))pcaDat = dat[,colnames(dat) %in% props]
pheno = dat[,!colnames(dat) %in% props]
h = as.matrix(pcaDat)
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
f2_5 = ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)
f2_5# MEs = pcaDat
# plot = pheatmap::pheatmap(cor(MEs), border_color = 'white', kmeans_k = 3,
# fontsize = 7, treeheight_col = 0)
# plotDimensions
## [1] 73 37
Samples
row1 = cowplot::plot_grid(plotlist = list(f2_1, f2_2, f2_3, f2_4), ncol = 2, align = 'hv',
labels = c('b', 'c', 'd', 'e'), label_fontface = 'plain')
row2 = cowplot::plot_grid(plotlist = list(f2_5),
labels = c('f'), label_fontface = 'plain')
cowplot::plot_grid(plotlist = list(row1, row2), nrow = 2,
rel_heights = c(1, 1.2), axis = 'tb') if(exportFigures) {
ggsave('../manuscript/figures/2_Panel.pdf', width = 7, height = 10, dpi = 1200)
ggsave('../manuscript/figures/2_Panel.png', width = 7, height = 10, dpi = 1200)
}dat = obj$tests_combined_adj
pcaDat = dat[,colnames(dat) %in% props]
pheno = dat[,!colnames(dat) %in% props]
# pca = prcomp(pcaDat, center = T, scale. = F)
pca = prcomp(pcaDat, center = T, scale. = T)
PCs = colnames(as.data.frame(pca$x[,1:6]))
p = cbind(as.data.frame(pheno), as.data.frame(pca$x[,1:6]))
p$Sex = factor(p$Sex, levels = c('Female', 'Male'))
p$SexCondition = paste0(p$Sex, '_', dat$Condition)
p$SexCondition = factor(p$SexCondition,
levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
obj$analyses$PCA_OutcomeVars_BatchAdj = peigs <- pca$sdev^2
e = data.frame(
PC = colnames(pca$x),
Proportion = eigs/sum(eigs),
Cumulative = cumsum(eigs)/sum(eigs))
e = e[1:5,] %>%
mutate(PC = factor(PC, levels = colnames(pca$x)))
f3_1 = ggplot(e, aes(x = PC, y = Proportion*100)) +
geom_line(aes(group = 1), color = 'grey') +
geom_point(size = 2, pch = 21, fill = obj$colors[2]) +
scale_y_continuous(limits = c(5, 25), breaks = seq(0, 25, 5)) +
labs(x = '', y = 'Variance explained (%)')
f3_1Check Normality
idVars = c(colnames(pheno), 'SexCondition')
PCs = c('PC1', 'PC2', 'PC3', 'PC4')
pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
dplyr::select(idVars, PCs) %>%
gather(PC, value, -idVars) %>%
group_by(PC, Sex, Condition) %>% nest() %>%
mutate(SW = map(data, function(d) shapiro.test(d$value)),
tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()Normality is not violated in any groups
Check heteroskedasticity
pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
dplyr::select(idVars, PCs) %>%
gather(PC, value, -idVars) %>%
group_by(PC) %>% nest() %>%
dplyr::mutate(
t = map(data, function(d) {
car::leveneTest(value ~ Sex * Condition, d)
})
) %>% dplyr::select(-data) %>% unnest()
ppAll variances are homogenous
idVars = c(colnames(pheno), 'SexCondition')
PCs = c('PC1', 'PC2', 'PC3', 'PC4')
pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
dplyr::select(idVars, PCs) %>%
gather(PC, value, -idVars) %>%
group_by(PC) %>% nest() %>%
mutate(
m = map2(data, PC, function(d, p) {
summary(aov(value ~ Condition * Sex, d))
})
) %>% dplyr::select(-data)
plots = obj$analyses$PCA_OutcomeVars_BatchAdj %>%
dplyr::select(idVars, PCs) %>%
gather(PC, value, -idVars) %>%
group_by(PC) %>% nest() %>% mutate(
plot = map2(data, PC, function(d, p) {
ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition),
position = position_dodge(0.6)) +
geom_point(pch = 21, position = position_dodge(0.6), size = 2.5) +
labs(x = '', y = p) + facet_grid(~Sex) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
theme(legend.position = 'none')
})
)PC1:
## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 1 3.3 3.304 0.411 0.523
## Sex 1 0.1 0.140 0.017 0.895
## Condition:Sex 1 8.4 8.401 1.046 0.310
## Residuals 69 554.3 8.033
PC2:
## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 1 3.34 3.34 1.233 0.271
## Sex 1 139.50 139.50 51.449 0.000000000646 ***
## Condition:Sex 1 6.04 6.04 2.226 0.140
## Residuals 69 187.09 2.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PC3:
## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 1 28.70 28.699 7.897 0.00644 **
## Sex 1 1.10 1.097 0.302 0.58449
## Condition:Sex 1 0.87 0.874 0.241 0.62530
## Residuals 69 250.74 3.634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PC4:
## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 1 0.42 0.417 0.159 0.69176
## Sex 1 27.24 27.237 10.344 0.00198 **
## Condition:Sex 1 30.08 30.076 11.423 0.00120 **
## Residuals 69 181.68 2.633
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID)))
pp = left_join(baseline, obj$analyses$PCA_OutcomeVars_BatchAdj)
pp$Condition = factor(pp$Condition, levels = c('Control', 'CMS'))
f3_3 = ggplot(pp, aes(x = DS, y = PC1)) +
geom_smooth(method = 'lm', se = F, lty = 2, aes(col = SexCondition)) +
geom_point(pch = 21, size = 2, aes(fill = SexCondition)) +
labs(y = "PC1 score", x = "Baseline David's Score") +
facet_grid(Condition~Sex) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
theme(legend.position = 'none')
f3_3Stats
##
## Call:
## lm(formula = PC1 ~ DS * Sex, data = pp[pp$Condition == "CMS",
## ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9194 -1.7251 -0.1317 1.6727 5.5421
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.421 2.045 -2.162 0.03754 *
## DS 2.856 1.287 2.220 0.03301 *
## SexMale 9.246 2.591 3.568 0.00107 **
## DS:SexMale -5.834 1.637 -3.563 0.00108 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.554 on 35 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.2857, Adjusted R-squared: 0.2245
## F-statistic: 4.667 on 3 and 35 DF, p-value: 0.007594
## Df Sum Sq Mean Sq F value Pr(>F)
## DS 1 6.22 6.22 0.954 0.33533
## Sex 1 2.30 2.30 0.353 0.55639
## DS:Sex 1 82.78 82.78 12.695 0.00108 **
## Residuals 35 228.22 6.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 9 observations deleted due to missingness
p = obj$analyses$PCA_OutcomeVars_BatchAdj
PC = 'PC1'
p$PC = p[[PC]]
cors = sapply(1:ncol(pcaDat), function(prop) { cor(p$PC, pcaDat[,prop], method = 'spearman') })
pvals = sapply(1:ncol(pcaDat), function(prop) { cor.test(p$PC, pcaDat[,prop], method = 'spearman')$p.value })
cors = data.frame(feature = colnames(pcaDat), rho = cors, pval = pvals)
cors$feature = factor(cors$feature, levels = cors$feature[order(abs(cors$rho), decreasing = F)])
cors = cors[order(abs(cors$rho), decreasing = T), ]
cors$Sign = ifelse(cors$rho < 0, yes = 'neg', no = 'pos')
cors$Sign = factor(cors$Sign, levels = c('pos', 'neg'))
# cors$Sig = cors$pval <= 0.05
# cors$Sig = p.adjust(cors$pval, method = 'BH') <= 0.1
# ? Bonferroni under 0.1 - is that good enough ?
cors$Sig = p.adjust(cors$pval, method = 'bonferroni') <= 0.1
sum(cors$Sig)## [1] 16
cors$feature_pretty = plyr::mapvalues(cors$feature, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
f3_4 = ggplot(cors, aes(y = abs(rho), x = feature_pretty, col = Sign, fill = Sign)) +
geom_bar(stat = 'identity', width = 0.01) +
geom_point(size = 2, pch = 21) +
geom_point(data = subset(cors, Sig), size = 2.5, pch = 21, color = 'black') +
scale_fill_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '') +
scale_color_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '') +
# xlab('') + ylab() + coord_flip() + obj$theme +
xlab('') + ylab("Correlation") + coord_flip() + obj$theme +
theme(plot.title = element_text(hjust = 0.5), legend.position = 'top')
f3_4dat = obj$tests_combined_adj
idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber',
'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID))) %>%
rename(Dominance = DS)
d = left_join(dat, baseline)
idVars = c(idVars, 'Dominance')
plots = d %>%
dplyr::select(idVars, props) %>%
gather(prop, value, -idVars) %>%
group_by(prop)
plots$prop = plyr::mapvalues(plots$prop, props, props_pretty)
plots$Sex = factor(plots$Sex, levels = c('Female', 'Male'))
plots$Condition = factor(plots$Condition, levels = c('Control', 'CMS'))
plots$SexCondition = paste0(plots$Sex, '_', dat$Condition)
plots$SexCondition = factor(plots$SexCondition,
levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
selected_props = c('OFT | Distance', 'EPM | Closed arm time')
# selected_props = c('OFT | Inner zone distance', 'EPM | Open arm distance')
res = lapply(selected_props, function(i) {
pp = subset(plots, prop == i)
l = c(min(pp$value), max(pp$value))
p1 = ggplot(pp, aes(x = Condition, y = value, fill = SexCondition)) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
geom_point(pch = 21, size = 2,
position = position_jitterdodge(jitter.width = 0,
jitter.height = 0, dodge.width = 0.5)) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
facet_grid(Sex~.) + labs(x = '', y = i) +
scale_y_continuous(limits = l) +
theme(legend.position = 'none') +
theme(strip.text.y = element_blank())
p2 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Dominance, y = value, fill = Sex, color = Sex)) +
geom_smooth(method = 'lm', se = F, lty = 2) +
geom_point(pch = 21, size = 2, color = 'grey20') +
scale_fill_manual(values = obj$four_color[3:4]) +
scale_color_manual(values = obj$four_color[3:4]) +
scale_y_continuous(limits = l) +
labs(y = '', x = "Baseline David's Score") +
facet_grid(Sex~Condition) +
theme(legend.position = 'none',
axis.text.y = element_blank(),
plot.margin = margin(0,0,0,-10, "pt")
)
cowplot::plot_grid(plotlist = list(p1, p2), axis = 'tb',
align = 'h', rel_widths = c(1,1.5))
})
f3_5 = cowplot::plot_grid(plotlist = res, ncol = 2, align = 'v', rel_widths = c(1,1),
labels = c('e', 'f'), label_fontface = 'plain')
f3_5row1 = cowplot::plot_grid(plotlist = list(f3_1, f3_2), ncol = 2, axis = 't', align = 'h',
labels = c('a', 'b'), label_fontface = 'plain', rel_widths = c(0.7,1))
row2 = cowplot::plot_grid(plotlist = list(f3_3), ncol = 1, align = 'hv',
labels = c('c'), label_fontface = 'plain')
col1 = cowplot::plot_grid(plotlist = list(row1, row2), ncol = 1, axis = 't', align = 'h',
rel_widths = c(0.7,1), rel_heights = c(1, 2.1))
col2 = cowplot::plot_grid(plotlist = list(f3_4), ncol = 1, align = 'hv',
labels = c('d'), label_fontface = 'plain')
l1 = cowplot::plot_grid(plotlist = list(col1, col2), align = 'hv', ncol = 2, rel_widths = c(1, 0.9),
label_fontface = 'plain')
cowplot::plot_grid(plotlist = list(l1, f3_5), align = 'hv', ncol = 1, rel_heights = c(1, 0.6))idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber',
'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID))) %>%
rename(Dominance = DS) %>% group_by(GroupNumber, Sex) %>%
mutate(
Rank = rank(-Dominance),
Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
)
baseline = left_join(baseline, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])
p <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
dplyr::select(RealMouseNumber, PC1)
d = left_join(baseline, p)
d$Sex = factor(d$Sex, levels = c('Female', 'Male'))
d$Condition = factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition = paste0(d$Sex, '_', dat$Condition)
d$SexCondition = factor(d$SexCondition,
levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
d$Rank_3 = plyr::mapvalues(d$Rank, c('alpha', 'beta', 'gamma', 'delta'),
c('alpha', 'beta / gamma', 'beta / gamma', 'delta'))
d$Rank_3 = factor(d$Rank_3, levels = c('alpha', 'beta / gamma', 'delta'))
d$Rank_3_ordinal = as.numeric(as.character(plyr::mapvalues(d$Rank_3,
c('alpha', 'beta / gamma', 'delta'),
c(1, 0, -1))))
pp = d
pp$value = pp$PC1
l = c(min(pp$value), max(pp$value))
xx = data.frame(
Sex = c('Female', 'Female', 'Male', 'Male'),
x = c(1.5, 2.5, 1.5, 2.5),
y = c(5, 0, 4.5, 5),
Sig = c('*', '', '#', '*'),
Line = c('_', '', '_', '_')
)
f4_1 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Rank_3, y = value, fill = Rank_3)) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = Rank_3)) +
geom_point(pch = 21, size = 2,
position = position_jitterdodge(jitter.width = 0,
jitter.height = 0, dodge.width = 0.5)) +
geom_text(data = xx[-3,], inherit.aes = F, aes(x = x, y = y, label = Sig), fontface = 'bold', size = 4,
color = 'grey30') +
geom_text(data = xx[3,], inherit.aes = F, aes(x = x, y = y, label = Sig), size = 3, fontface = 'bold', color = 'grey30') +
geom_text(data = xx, inherit.aes = F, aes(x = x, y = y, label = Line), fontface = 'bold', color = 'grey30') +
scale_color_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
scale_fill_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4]),
labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
scale_x_discrete(labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
scale_y_continuous(limits = l) +
facet_grid(Sex~Condition) + labs(x = '', y = 'PC1 score') +
theme(legend.position = 'none',
axis.title.x = element_text(family="Helvetica")
)
f4_1dat = subset(pp, Condition == 'CMS')
dat$value = dat[['PC1']]
dat %>% group_by(Sex, Rank_3) %>%
filter(!is.na(value)) %>%
dplyr::summarize(
n = n(),
n_groups = length(unique(GroupNumber)),
M = mean(value, na.rm = T),
sd = sd(value, na.rm = T)
)dat %>% group_by(Sex, Rank_3) %>% nest() %>%
mutate(SW = map(data, function(d) shapiro.test(d$value)),
tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()Normality is not violated
| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 5 | 0.8531606 | 0.5225123 |
| 33 | NA | NA |
Groupwise variances are homogenous
Test:
Males - omnibus - rank ordinal
Alpha to beta/gamma
dds = dat %>%
filter(Sex == 'Male' & Rank_3_ordinal %in% c(0, 1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])Delta to beta/gamma
dds = dat %>%
filter(Sex == 'Male' & Rank_3_ordinal %in% c(0, -1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])Females - omnibus
Alpha to beta/gamma
dds = dat %>%
filter(Sex == 'Female' & Rank_3_ordinal %in% c(0, 1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])Delta to beta/gamma
dds = dat %>%
filter(Sex == 'Female' & Rank_3_ordinal %in% c(0, -1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])Select top PC1 behaviors
# selected_props = as.character(cors$feature[abs(cors$rho) >= 0.5])
selected_props = as.character(cors$feature[cors$Sig == T])
dat = obj$tests_combined_adj
props = colnames(dat)[regexpr('_', colnames(dat)) > 0]
pcaDat = dat[,colnames(dat) %in% props]
pheno = dat[,!colnames(dat) %in% props]
pcaDat = pcaDat[,colnames(pcaDat) %in% selected_props]
h = as.matrix(pcaDat)
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)All Controls
h = as.matrix(pcaDat[pheno$Condition == 'Control',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)All CMS
h = as.matrix(pcaDat[pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)h = as.matrix(pcaDat[pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
# plot(hc, hang = -1, cex = 0.6)
hc$labels = plyr::mapvalues(hc$labels, props, props_pretty)
hcd <- as.dendrogram(hc)
# Default plot
# plot(hcd, type = "rectangle", ylab = "Height")
ggdendro::ggdendrogram(hc, rotate = T)#
# dend <- as.dendrogram(hc)
# # Extract the data (for rectangular lines)
# # Type can be "rectangle" or "triangle"
# dend_data <- ggdendro::dendro_data(dend, type = "rectangle")
# ggplot(dend_data$segments) +
# geom_segment(aes(x = x, y = y, xend = xend, yend = yend))Male Controls
h = as.matrix(pcaDat[pheno$Sex == 'Male' & pheno$Condition == 'Control',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)Male CMS
h = as.matrix(pcaDat[pheno$Sex == 'Male' & pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)Female Control
h = as.matrix(pcaDat[pheno$Sex == 'Female' & pheno$Condition == 'Control',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)Female CMS
h = as.matrix(pcaDat[pheno$Sex == 'Female' & pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = 'white') +
scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
scale_y_discrete(position = 'right') +
labs(x = '', y = '') +
coord_fixed() +
theme(
axis.text.x = element_blank(), panel.border = element_blank(),
plot.margin = margin(15,0,0,15, "pt"),
legend.position = 'bottom'
)dat = obj$tests_combined_adj
idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber',
'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID))) %>%
rename(Dominance = DS) %>% group_by(GroupNumber, Sex) %>%
mutate(
Rank = rank(-Dominance),
Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
)
d = left_join(dat, baseline)
idVars = c(idVars, 'Dominance', 'Rank')
d = d %>%
dplyr::select(idVars, props) %>%
gather(prop, value, -idVars) %>%
group_by(prop)
d$prop = plyr::mapvalues(d$prop, props, props_pretty)
d$Sex = factor(d$Sex, levels = c('Female', 'Male'))
d$Condition = factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition = paste0(d$Sex, '_', dat$Condition)
d$SexCondition = factor(d$SexCondition,
levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
selected_props = as.character(cors$feature_pretty[cors$Sig == T])
d = d %>% filter(prop %in% selected_props)
d$Rank_3 = plyr::mapvalues(d$Rank, c('alpha', 'beta', 'gamma', 'delta'),
c('alpha', 'beta / gamma', 'beta / gamma', 'delta'))
d$Rank_3 = factor(d$Rank_3, levels = c('alpha', 'beta / gamma', 'delta'))
d$Rank_3_ordinal = as.numeric(as.character(plyr::mapvalues(d$Rank_3,
c('alpha', 'beta / gamma', 'delta'),
c(1, 0, -1))))r1 = d %>% group_by(prop, Sex) %>%
nest() %>%
mutate(
res = map(data, function(dds) {
dds = dds[!(dds$Condition == 'CMS' & dds$Rank_3_ordinal %in% c(0, -1)),]
dds = dds %>%
ungroup() %>%
dplyr::mutate(
value_ranked = rank(value)
)
# x = as.data.frame(summary(lmPerm::aovp(value_ranked ~ Condition, dds))[[1]])
x = as.data.frame(summary(lmPerm::aovp(value ~ Condition, dds))[[1]])
x = x[1,]
# x = tidy(kruskal.test(value ~ Condition, dds))
x$comparison = 'alpha_to_control'
return(x)
})
) %>% dplyr::select(-data) %>% unnest()## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
r2 = d %>% group_by(prop, Sex) %>%
nest() %>%
mutate(
res = map(data, function(dds) {
dds = dds[!(dds$Condition == 'CMS' & dds$Rank_3_ordinal %in% c(0, 1)),]
dds = dds %>%
ungroup() %>%
dplyr::mutate(
value_ranked = rank(value)
)
# x = as.data.frame(summary(lmPerm::aovp(value_ranked ~ Condition, dds))[[1]])
x = as.data.frame(summary(lmPerm::aovp(value ~ Condition, dds))[[1]])
x = x[1,]
# x = tidy(kruskal.test(value ~ Condition, dds))
x$comparison = 'delta_to_control'
return(x)
})
) %>% dplyr::select(-data) %>% unnest()## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
r %>% filter(comparison == 'alpha_to_control' & Sex == 'Male') %>%
dplyr::select(comparison, Sex, prop, p.value)%>%
mutate(padj = p.adjust(p.value, method = 'BH')) %>% arrange(p.value)r %>% filter(comparison == 'alpha_to_control' & Sex == 'Female') %>%
dplyr::select(comparison, Sex, prop, p.value)%>%
mutate(padj = p.adjust(p.value, method = 'BH')) %>% arrange(p.value)r %>% filter(comparison == 'delta_to_control' & Sex == 'Male') %>%
dplyr::select(comparison, Sex, prop, p.value)%>%
mutate(padj = p.adjust(p.value, method = 'BH')) %>% arrange(p.value)r %>% filter(comparison == 'delta_to_control' & Sex == 'Female') %>%
dplyr::select(comparison, Sex, prop, p.value)%>%
mutate(padj = p.adjust(p.value, method = 'BH')) %>% arrange(p.value)rr = r %>%
filter(Sex == 'Male' & comparison == 'alpha_to_control') %>%
arrange(-p.value) %>%
dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr1 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
coord_flip() + theme(legend.position = 'none') +
ggtitle('CMS \u03B1-males vs. Controls')
rr = r %>%
filter(Sex == 'Male' & comparison == 'delta_to_control') %>%
arrange(-p.value) %>%
dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr2 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
coord_flip() + theme(legend.position = 'none') +
ggtitle('CMS \u03B4-males vs. Controls')
rr = r %>%
filter(Sex == 'Female' & comparison == 'alpha_to_control') %>%
arrange(-p.value) %>%
dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr3 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
coord_flip() + theme(legend.position = 'none') +
ggtitle('CMS \u03B1-females vs. Controls')
cowplot::plot_grid(plotlist = list(rr1, rr2, rr3), ncol = 3)extra:
rr = r %>%
filter(Sex == 'Female' & comparison == 'delta_to_control') %>%
arrange(-p.value) %>%
dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
ggplot(rr, aes(x = prop, y = -log10(p.value))) +
geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
coord_flip() + theme(legend.position = 'none') +
ggtitle('CMS \u03B1-females vs. Controls')rr = r %>% filter(prop %in% c('OFT | Distance', 'EPM | Closed arm time')) %>%
# rr = r %>% filter(prop %in% c('EPM | Open arm distance', 'EPM | Closed arm time')) %>%
dplyr::mutate(Sig = p.value <=.05)
rr$Sig = ifelse(rr$p.value >.1, NA,
ifelse(rr$p.value >= 0.05, '#',
ifelse(rr$p.value >= .01, '*',
ifelse(rr$p.value >= 0.001, "**", '***'))))
res = d %>% group_by(prop) %>% nest() %>%
filter(prop %in% c('OFT | Distance', 'EPM | Closed arm time')) %>%
mutate(
plot = map2(data, prop, function(pp, i) {
l = c(min(pp$value), max(pp$value))
p1 = ggplot(pp[pp$Condition == 'Control',], aes(x = as.factor(1), y = value, fill = Rank_3)) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, fill = 'grey', color = 'grey') +
geom_point(pch = 21, size = 2, position = position_jitter(height = 0, width = 0)) +
scale_y_continuous(limits = l) +
scale_color_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
scale_fill_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
facet_grid(Sex~Condition) + labs(x = '', y = i) +
theme(legend.position = 'none',
strip.text.y = element_blank(),
axis.text.x = element_blank()
)
xx = subset(rr, prop == i)
xx$y = l[2]
xx$Rank_3 = plyr::mapvalues(xx$comparison, c('alpha_to_control', 'delta_to_control'),
c('alpha', 'delta'))
p2 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Rank_3, y = value, fill = Rank_3)) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = Rank_3)) +
geom_point(pch = 21, size = 2,
position = position_jitterdodge(jitter.width = 0,
jitter.height = 0, dodge.width = 0.5)) +
scale_color_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
scale_fill_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4]),
labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
scale_x_discrete(labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
scale_y_continuous(limits = l) +
geom_text(data = xx, inherit.aes = F, aes(x = Rank_3, y = y, label = Sig), fontface = 'bold', color = 'grey30') +
facet_grid(Sex~Condition) + labs(x = '', y = '') +
theme(legend.position = 'none',
plot.margin = margin(0,0,0,-18, "pt"),
axis.text.y = element_blank(),
axis.title.x = element_text(family="Helvetica")
)
cowplot::plot_grid(plotlist = list(p1, p2), axis = 'tb',
align = 'h', rel_widths = c(1, 1.3))
})
)cowplot::plot_grid(plotlist = list(f4_1, res$plot[[2]], res$plot[[1]]), ncol = 3, align = 'hv',
rel_widths = c(0.8, 1, 1),
labels = c('a', 'b', 'c'), label_fontface = 'plain')wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = wins, nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F)
h = h %>% dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID)))
h = left_join(h, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])
p = subset(obj$profiles, !Suspect & Day %in% 2:4) %>%
dplyr::select(RealMouseNumber, all_props) %>%
group_by(RealMouseNumber) %>%
dplyr::summarise_at(.vars = all_props, mean, na.rm = T)
d = left_join(h, p)
idVars = c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props = all_props
p <- d %>%
dplyr::select(idVars, props) %>%
gather(prop, value, -idVars) %>%
group_by(prop, Sex) %>% nest() %>%
mutate(
M = map(data, function(d) {
tidy(cor.test(d$value, d$DS, method = 'spearman'))
})
) %>% dplyr::select(-data) %>% unnest() %>%
mutate(
log10pval = -log10(p.value),
prop = fct_reorder(prop, abs(estimate)),
sign = ifelse(estimate > 0, '+', '-'),
sig_nominal = p.value <= 0.05
) %>% group_by(Sex) %>%
mutate(qval = p.adjust(p.value, method = 'BH')) %>%
arrange(qval)
plots = d %>%
dplyr::select(idVars, props) %>%
gather(prop, value, -idVars) %>%
group_by(prop) %>% nest() %>%
mutate(
plot = pmap(list(data, prop),
function(d, p) {
ggplot(d, aes(x = DS, y = value, fill = Sex)) +
geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
geom_point(pch = 21, size = 2) +
labs(y = p, x = "Baseline David's Score") +
facet_grid(.~Sex) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
theme(legend.position = 'none')
})) %>% dplyr::select(-data)
p = left_join(p, plots)p = subset(p, !is.na(estimate))
p$rho = abs(p$estimate)
p$rho[p$qval > 0.1] = NA
p = subset(p, qval <= 0.1)
p$sign = factor(p$sign, levels = c('+', '-'))
p$prop_pretty = plyr::mapvalues(p$prop, all_props, all_props_pretty)
ord = p$prop_pretty[p$Sex == 'Male'][order(p$estimate[p$Sex == 'Male'], decreasing = F)]
p$prop_pretty = factor(p$prop_pretty, levels = ord)
p = subset(p, !is.na(prop_pretty))
p$Sex = factor(p$Sex, levels = c('Male', 'Female'))
hm.palette <- colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(p[!is.na(p$sign),], aes(y = Sex, x = prop_pretty, col = sign, fill = sign)) +
geom_point(pch = 22, aes(size = rho)) +
scale_fill_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = 'Direction') +
scale_color_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = 'Direction') +
scale_size_continuous(name = 'Correlation') +
labs(x = '', y = '') + coord_flip() +
theme(plot.title = element_text(hjust = 0.5), legend.position = 'right',
panel.border = element_blank())Hierarchy properties based on chases
Daily & mean of days 1:4 and day 5 separately
domProps = c('Steepness', 'Despotism', 'DirectionalConsistency', 'Landau')
domPropsPretty = c('Steepness', 'Despotism', 'Directional Consistency', "Landau's modified h'")
wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:5), wins = wins, nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = T)
h$uGroup = paste0(h$GroupNumber, h$Sex)
plots = lapply(1:(length(domProps)-1), function(i) {
prop = domProps[i]
h$prop = h[[prop]]
lims = c(min(h$prop, na.rm = T), max(h$prop, na.rm = T))
p1 = ggplot(h, aes(x = as.factor(Day), y = prop, fill = Sex)) +
geom_line(aes(color = Sex, group = uGroup), lty = 1, size = 0.2, position = position_dodge(0)) +
geom_point(aes(group = uGroup, col = Sex), pch = 19, size = 1.2, position = position_dodge(0)) +
labs(x = '', y = domPropsPretty[i]) +
scale_color_manual(values = obj$colors_sex) +
scale_fill_manual(values = obj$colors_sex) +
scale_y_continuous(limits = lims) +
theme(legend.position = 'none', axis.title.x=element_blank())
p2 = h %>% dplyr::filter(Day %in% 1:4) %>%
group_by(GroupNumber, Condition, Sex) %>%
dplyr::summarise(prop = mean(prop, na.rm = T)) %>%
ggplot(aes(x = Sex, y = prop, fill = Sex)) +
geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
labs(x = '', y = '') +
scale_y_continuous(limits = lims) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
theme(legend.position = 'none',
axis.text.y = element_blank(),
axis.title=element_blank()
) +
ggtitle('Baseline')
p3 = h %>% dplyr::filter(Day %in% 5) %>%
group_by(GroupNumber, Condition, Sex) %>%
dplyr::summarise(prop = mean(prop, na.rm = T)) %>%
ggplot(aes(x = Sex, y = prop, fill = Sex)) +
geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
labs(x = '', y = '') +
scale_y_continuous(limits = lims) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
theme(legend.position = 'none',
axis.text.y = element_blank(),
axis.title=element_blank()
) +
ggtitle('Acute stress')
cowplot::plot_grid(plotlist = list(p1, p2, p3), ncol = 3, align = 'h',
label_fontface = 'plain',
rel_widths = c(3.4, 1, 1))
})
lastPlot = lapply(length(domProps), function(i) {
prop = domProps[i]
h$prop = h[[prop]]
lims = c(min(h$prop, na.rm = T), max(h$prop, na.rm = T))
p1 = ggplot(h, aes(x = as.factor(Day), y = prop, fill = Sex)) +
geom_boxplot(width = 0.7, alpha = 0.8, outlier.shape = NA, aes(color = Sex)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
labs(x = 'Day', y = domPropsPretty[i]) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
scale_y_continuous(limits = lims) +
theme(legend.position = 'none')
p2 = h %>% dplyr::filter(Day %in% 1:4) %>%
group_by(GroupNumber, Condition, Sex) %>%
dplyr::summarise(prop = mean(prop, na.rm = T)) %>%
ggplot(aes(x = Sex, y = prop, fill = Sex)) +
geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
# geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.3) +
labs(x = '', y = '') +
scale_y_continuous(limits = lims) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
theme(legend.position = 'none', axis.text.y = element_blank(),
axis.title.y=element_blank()) +
ggtitle('Baseline')
p3 = h %>% dplyr::filter(Day %in% 5) %>%
group_by(GroupNumber, Condition, Sex) %>%
dplyr::summarise(prop = mean(prop, na.rm = T)) %>%
ggplot(aes(x = Sex, y = prop, fill = Sex)) +
geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.3) +
labs(x = '', y = '') +
scale_y_continuous(limits = lims) +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
theme(legend.position = 'none', axis.text.y = element_blank(),
axis.title.y=element_blank()) +
ggtitle('Acute stress')
cowplot::plot_grid(plotlist = list(p1, p2, p3), ncol = 3, align = 'h',
rel_widths = c(3.4, 1, 1))
})
plots[[length(domProps)]] = lastPlot[[1]]
cowplot::plot_grid(plotlist = plots, ncol = 1, align = 'v', axis = 'tblr',
labels = c('a', 'b', 'c', 'd'), label_fontface = 'plain',
rel_heights = c(rep(1, length(domProps)-1), 1.1))if(exportFigures) {
ggsave('../manuscript/figures/Supp_HierarchyPanel.pdf', width = 8.5, height = 11, dpi = 1200)
ggsave('../manuscript/figures/Supp_HierarchyPanel.png', width = 8.5, height = 11, dpi = 1200)
}Stats:
domProps = c('Steepness', 'Despotism', 'DirectionalConsistency', 'Landau')
domPropsPretty = c('Steepness', 'Despotism', 'Directional Consistency', "Landau's modified h'")
wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:5), wins = wins, nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = T)
h$uGroup = paste0(h$GroupNumber, h$Sex)
h$Stage = plyr::mapvalues(h$Day, 1:5, c(rep('Baseline', 4), 'Acute stress'))
h$Stage = factor(h$Stage, levels = c('Baseline', 'Acute stress'))
h = h %>% group_by(uGroup, GroupNumber, Sex, Stage) %>%
dplyr::summarise_at(.vars = domProps, mean, na.rm = T)
h %>% gather(prop, value, -uGroup, -GroupNumber, -Sex, -Stage) %>%
group_by(prop, Sex, Stage) %>% nest() %>%
dplyr::mutate(
SW = map(data, function(d) shapiro.test(d$value)),
tSW = map(SW, broom::tidy)
) %>%
select(-data, -SW) %>% unnest()Normality violated only for Landau’s h in female groups during acute stress
Check heteroskedasticity
h %>% gather(prop, value, -uGroup, -GroupNumber, -Sex, -Stage) %>%
group_by(prop) %>% nest() %>%
dplyr::mutate(
t = map(data, function(d) {
car::leveneTest(value ~ Sex * Stage, d)
})
)%>%
select(-data) %>% unnest()Variances are heteroskedastic for despotism
Test
##
## Error: uGroup
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.4634 0.4634 21.181 0.000194 ***
## Stage 1 0.0247 0.0247 1.131 0.300839
## Residuals 19 0.4157 0.0219
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 0.00001 0.000015 0.002 0.964
## Sex:Stage 1 0.00109 0.001091 0.157 0.696
## Residuals 19 0.13186 0.006940
##
## Error: uGroup
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.3859 0.3859 7.834 0.0111 *
## Residuals 20 0.9852 0.0493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 0.00707 0.007074 0.757 0.395
## Sex:Stage 1 0.00600 0.006003 0.642 0.432
## Residuals 20 0.18691 0.009345
##
## Error: uGroup
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.5689 0.5689 24.405 0.0000909 ***
## Stage 1 0.1325 0.1325 5.685 0.0277 *
## Residuals 19 0.4429 0.0233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 0.00723 0.00723 0.611 0.444
## Sex:Stage 1 0.00270 0.00270 0.228 0.638
## Residuals 19 0.22495 0.01184
Interesting effect of stage
Males:
##
## Error: uGroup
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 0.1325 0.13252 3.153 0.114
## Residuals 8 0.3362 0.04202
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 0.00901 0.009015 0.523 0.49
## Residuals 8 0.13779 0.017223
Females:
##
## Error: uGroup
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11 0.1067 0.009702
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 0.00092 0.000916 0.116 0.74
## Residuals 11 0.08716 0.007924
ggplot(h, aes(x = Stage, y = DirectionalConsistency, fill = Sex)) +
geom_line(aes(group = uGroup, col = Sex)) +
geom_point(size = 2, pch = 21) +
labs(x = '') +
scale_fill_manual(values = obj$colors_sex) +
scale_color_manual(values = obj$colors_sex) +
facet_grid(~Sex) +
theme(legend.position = 'none', axis.text.y = element_blank(),
axis.title.y=element_blank())##
## Error: uGroup
## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.1363 0.13632 2.113 0.162
## Stage 1 0.0709 0.07092 1.099 0.308
## Residuals 19 1.2260 0.06452
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Stage 1 0.0012 0.00121 0.050 0.825
## Sex:Stage 1 0.0384 0.03841 1.598 0.222
## Residuals 19 0.4568 0.02404
Baseline non-parametric:
##
## Kruskal-Wallis rank sum test
##
## data: Landau by Sex
## Kruskal-Wallis chi-squared = 5.8385, df = 1, p-value = 0.01568
Acute stress non-parametric:
##
## Kruskal-Wallis rank sum test
##
## data: Landau by Sex
## Kruskal-Wallis chi-squared = 0.022159, df = 1, p-value = 0.8817
idVars = c(colnames(pheno), 'SexCondition')
PCs = c('PC2', 'PC3')
# PCs = c('PC2', 'PC3', 'PC4', 'PC5')
pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
dplyr::select(idVars, PCs) %>%
gather(PC, value, -idVars) %>%
group_by(PC) %>% nest() %>%
mutate(
m = map2(data, PC, function(d, p) {
summary(aov(value ~ Condition * Sex, d))
})
) %>% dplyr::select(-data)
plots = obj$analyses$PCA_OutcomeVars_BatchAdj %>%
dplyr::select(idVars, PCs) %>%
gather(PC, value, -idVars) %>%
group_by(PC) %>% nest() %>% mutate(
plot = map2(data, PC, function(d, p) {
ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition),
position = position_dodge(0.6)) +
geom_point(pch = 21, position = position_dodge(0.6), size = 2.5) +
labs(x = '', y = p) + facet_grid(~Sex) +
scale_fill_manual(values = obj$four_color) +
scale_color_manual(values = obj$four_color) +
theme(legend.position = 'none')
})
)
pp$m[[1]]## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 1 3.34 3.34 1.233 0.271
## Sex 1 139.50 139.50 51.449 0.000000000646 ***
## Condition:Sex 1 6.04 6.04 2.226 0.140
## Residuals 69 187.09 2.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dat = obj$tests_combined_adj
idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber',
'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>%
dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>%
mutate(MouseID = as.numeric(as.character(MouseID))) %>%
rename(Dominance = DS) %>% group_by(GroupNumber, Sex) %>%
mutate(
Rank = rank(-Dominance),
Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
)
d = left_join(dat, baseline)
idVars = c(idVars, 'Dominance', 'Rank')
d = d %>%
dplyr::select(idVars, props) %>%
gather(prop, value, -idVars) %>%
group_by(prop)
d$prop = plyr::mapvalues(d$prop, props, props_pretty)
d$Sex = factor(d$Sex, levels = c('Female', 'Male'))
d$Condition = factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition = paste0(d$Sex, '_', dat$Condition)
d$SexCondition = factor(d$SexCondition,
levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
selected_props = as.character(cors$feature_pretty[cors$Sig == T])
d = d %>% filter(prop %in% selected_props)
d$Rank_3 = plyr::mapvalues(d$Rank, c('alpha', 'beta', 'gamma', 'delta'),
c('alpha', 'beta / gamma', 'beta / gamma', 'delta'))
d$Rank_3 = factor(d$Rank_3, levels = c('alpha', 'beta / gamma', 'delta'))
d$Rank_3_ordinal = as.numeric(as.character(plyr::mapvalues(d$Rank_3,
c('alpha', 'beta / gamma', 'delta'),
c(1, 0, -1))))
r1 = d %>% group_by(prop, Sex) %>%
nest() %>%
mutate(
res = map(data, function(dds) {
dds = dds[!(dds$Condition == 'CMS' & dds$Rank_3_ordinal %in% c(0, -1)),]
# dds = dds %>%
# ungroup() %>%
# dplyr::mutate(
# value_ranked = rank(value)
# )
# # x = as.data.frame(summary(lmPerm::aovp(value_ranked ~ Condition, dds))[[1]])
x = as.data.frame(summary(lmPerm::aovp(value ~ Condition, dds))[[1]])
x = x[1,]
xx = dds %>% group_by(Condition) %>%
dplyr::summarise(M = mean(value, na.rm = T))
x$direction = ifelse(xx$M[xx$Condition == 'CMS'] > xx$M[xx$Condition == 'Control'], 'positive', 'negative')
# x = tidy(kruskal.test(value ~ Condition, dds))
x$comparison = 'alpha_to_control'
return(x)
})
) %>% dplyr::select(-data) %>% unnest()## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
r2 = d %>% group_by(prop, Sex) %>%
nest() %>%
mutate(
res = map(data, function(dds) {
dds = dds[!(dds$Condition == 'CMS' & dds$Rank_3_ordinal %in% c(0, 1)),]
dds = dds %>%
ungroup() %>%
dplyr::mutate(
value_ranked = rank(value)
)
# x = as.data.frame(summary(lmPerm::aovp(value_ranked ~ Condition, dds))[[1]])
x = as.data.frame(summary(lmPerm::aovp(value ~ Condition, dds))[[1]])
x = x[1,]
xx = dds %>% group_by(Condition) %>%
dplyr::summarise(M = mean(value, na.rm = T))
x$direction = ifelse(xx$M[xx$Condition == 'CMS'] > xx$M[xx$Condition == 'Control'], 'positive', 'negative')
# x = tidy(kruskal.test(value ~ Condition, dds))
x$comparison = 'delta_to_control'
return(x)
})
) %>% dplyr::select(-data) %>% unnest()## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
## [1] "Settings: unique SS "
rr = r %>%
filter(Sex == 'Male' & comparison == 'alpha_to_control') %>%
arrange(-p.value) %>%
dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr1 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
coord_flip() + theme(legend.position = 'none') +
ggtitle('CMS \u03B1-males vs. Controls')
rr = r %>%
filter(Sex == 'Male' & comparison == 'delta_to_control') %>%
arrange(-p.value) %>%
dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr2 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
coord_flip() + theme(legend.position = 'none') +
ggtitle('CMS \u03B4-males vs. Controls')
rr = r %>%
filter(Sex == 'Female' & comparison == 'alpha_to_control') %>%
arrange(-p.value) %>%
dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr3 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
coord_flip() + theme(legend.position = 'none') +
ggtitle('CMS \u03B1-females vs. Controls')
cowplot::plot_grid(plotlist = list(rr1, rr2, rr3), ncol = 3)